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1  Introduction 

The  research  funded  under  this  grant  focused  on  solution  algorithms  for 
porous  media  flow  models.  We  sought  to  improve  both  spatial  and  temporal 
computational  methods  for  a  range  of  problems  associated  with  groundwater 
and  surface  water  flows. 

Groundwater  flow  problems  can  be  difficult  to  resolve  for  several  reasons. 
The  governing  equations  are  nonlinear,  as  the  permeability  coefficient  for 
the  Darcy  velocity  depends  on  the  saturation  level.  This  dependence  on 
saturation  also  means  that  the  governing  equation  can  change  type,  from 
parabolic  to  hyperbolic,  or  vice  versa,  as  an  infiltration  front  moves  through 
the  domain.  In  addition,  the  governing  equation  can  change  type  as  fluid 
flows  from  one  material  type  to  another  (e.g.,  from  sand  to  clay).  In  all 
cases,  the  challenges  are  both  on  the  temporal  and  spatial  domains. 

2  Unsaturated  Flow 

Richards’  Equation  is  still  one  of  the  dominant  modeling  equations  for  va- 
dose  zone  flow.  The  equation  can  be  formulated  with  either  water  pressure 
or  water  saturation  as  the  primary  unknown,  or  the  equation  can  be  posed 
in  mixed  form  where  both  variables  appear  as  unknowns.  While  modeling 
vadose  zone  behavior  is  increasingly  simulated  using  two-phase  flow  equa¬ 
tions,  resolution  of  Richards’  Equation  remains  popular  because  the  size  of 
the  systems  is  significantly  smaller  than  that  for  the  two-phase  equations. 
In  addition,  the  two  phase  equations  include  model  equations  for  a  com¬ 
pressible  fluid  (i.e.,  air),  which  pose  computational  difficulties  as  the  density 
of  the  fluid  now  depends  on  the  air  pressure. 

Accurate  solution  of  Richards’  Equation  requires  resolution  in  both  space 
and  time.  Both  spatial  and  temporal  adaptive  routines  have  been  studied  for 
this  problem,  as  have  higher  order  discretization  methods  in  both  space  and 
time.  The  state-of-the-art  in  adaptivity  was  pushed  forward  for  Richards’ 
Equation  in  [11],  In  that  article,  Farthing  and  his  co-authors  designed  a 
space-time  error  measure  and  algorithm  for  Richards’  Equation.  They  base 
their  algorithm  on  a  loose  coupling  of  the  two  distinct  discretizations  to 
enable  users  to  take  advantage  of  existing,  mature  spatial  and  temporal 
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adaption  technology.  This  approach  is  important  for  ARO  applications  as 
extensive  effort  has  been  invested  in  developing  flow  simulation  codes.  The 
loose  coupling  allows  these  codes  to  be  reused  with  significant  improvements 
in  accuracy  and  efficiency. 

The  space/tinre  adaptive  algorithm  was  extended  in  [9],  [10]  to  use  local 
discontinuous  Galerkin  methods  for  spatial  discretization.  The  discontin¬ 
uous  Galerkin  method  is  attractive  because  of  it  ability  to  conserve  mass 
locally,  to  increase  the  order  of  the  approximation  in  specific  regions  of  the 
domain,  to  accurately  approximate  solutions  that  are  discontinuous  or  nearly 
discontinuous,  and  to  parallelize  easily.  The  primary  drawback  is  that  there 
are  a  large  number  of  degrees  of  freedom  associated  with  the  discrete  system. 
This  difficulty  is  overcome  by  using  discontinuous  Galerkin  methods  only  in 
subregions  of  the  domain  that  benefit  from  their  use,  which  is  described  and 
analyzed  in  [9],  [10]. 

3  Heterogeneous  Media 

We  developed  a  test  set  of  real-word,  physical  application  problems  that 
would  exhibit  the  numerical  difficulties  associated  with  their  simulation. 
The  set  included  problems  from  the  literature  and  problems  that  are  inher¬ 
ently  ill-conditioned,  yet  scientifically  important.  The  document  included 
descriptions  of: 

•  infiltration  problems  with  increases  in  levels  of  the  water  table; 

•  trenches  lined  with  relatively  impermeable  materials,  much  like  land¬ 
fills; 

•  flow  throughly  media  with  variably  distributed  heterogeneity. 

Much  of  the  work  centered  around  resolution  of  Richards’  Equation,  but 
the  techniques  developed  to  handle  problems  in  this  test  set  can  be  extended 
to  handle  two-phase  flow  equations.  The  test  set  was  uploaded  and  described 
in  an  interim  progress  report. 

4  Coupled  Free  Flow  and  Darcy  Flow 

A  significant  amount  of  work  has  been  done  in  the  last  decade  to  develop 
algorithms  for  resolving  coupled  Stokes  and  Darcy  flows  (see,  for  instance, 
[8]  and  references  therein).  However,  this  work  has  involved  Newtonian  flu¬ 
ids,  where  the  extra  stress  tensor  is  directly  proportional  to  the  deformation 
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tensor.  The  constant  of  proportionality  is  the  dynamic  viscosity  of  the  fluid. 

We  have  examined  this  coupling  for  non-Newtonian  fluids,  where  the  vis¬ 
cosity  of  the  fluid  is  non-constant  and  depends  on  the  magnitude  of  the 
deformation  tensor  [7] 

5  Resolving  Coupled  Parabolic/Hyperbolic  Flow 

We  have  also  developed  algorithms  for  simple  post-processing  of  discrete, 
continuous  Galerkin  approximations  of  the  convection-diffusion  equations 
where  the  equation  type  changes  in  both  space  and  time.  Sharp  fronts 
appear  in  the  transition  region  as  the  equation  moves  from  hyperbolic  to 
parabolic,  and  numerical  approximations  along  the  front  exhibit  local,  non¬ 
physical  oscillations.  Upwinding  methods  have  been  used  to  resolve  this 
problem,  as  have  discontinuous  Galerkin  and  local  discontinuous  Galerkin 
methods.  Upwinding  methods  often  perform  poorly  when  the  grid  is  not 
aligned  with  the  primary  flow  direction,  and,  in  temporally  dependent  prob¬ 
lems,  the  operator  must  be  reconstructed  at  each  time  step.  As  mentioned 
earlier,  discontinous  Galerkin  methods  use  large  degrees  of  freedom,  and 
local  discontinuous  Galerkin  methods  require  remeshing  in  temporally  de¬ 
pendent  problems  where  the  sharp  front  moves  as  a  function  of  time  (e.g., 
infiltration  problems).  In  [6] 

In  this  example,  the  domain  is  parabolic  on  [0, 1),  hyperbolic  on  [1.5,2), 
and  parabolic  again  on  [2,2.5].  We  define  t  E  [0,2],  with  At  =  ygg.  The 
true  solution  is  given  by  the  function 

{—  sin(x  —  1)  +  exp(x  —  1)  +  t(x  —  l)2,  0  <  x  <  1 

(x  —  l)2(t  +  1)  —  1,  1  <  x  <  1.5 

1  +  (x  -  1.5)  (x2  +  (0.25*  -  3.5)  x  +  (2.75  -  0.625*))  ,  1.5  <  x  <  2 

The  solution  is  discontinuous  at  1.5,  but  the  flux  is  continuous  on  [0,2], 
The  unstabilized  solution  exhibits  highly  oscillatory  behavior,  but  both  of 
the  stabilized  solutions  remove  these  oscillations. 

The  approximation  is  computed  on  a  grid  with  h  =  p  and  At  =  j^g. 

In  the  two  grid  solution,  shown  on  the  right  in  Figure  2,  the  smoothing 
operator  is  computed  on  a  grid  with  h  =  . 
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Figure  1:  True  solution  (left),  unstabilized  solution  (right) 


1  1.2  1.4  1.6  1.8  2 


Figure  2:  Smoothed  solution  on  one  grid  (left)  and  two  grids  (right) 


We  have  also  applied  the  operator  to  one-  and  two-dimensional  infiltra¬ 
tion  problems.  The  one-dimensional  results  are  shown  in  Figures  3  and  4, 
and  the  two-dimensional  results  are  in  [6].  The  one-dimensional  problem 
shows  the  solution  to  the  governing  equation 

ut  +  uux  —  0.001  uxx  =  /,  —  1  <  x  <  4,  t  >  0. 

where  /  is  defined  so  that  the  true  solution  is 

u(x,  t)  =  1  —  tanh  (  - (x  —  t) 

v  ’  VO-002  v  ’ 

The  true  solution  is  shown  on  the  left  in  Figure  3,  and  the  unstabilized 
solution  is  shown  on  the  right. 

The  results  on  the  left  in  Figure  4  were  computed  by  solving  the  convection- 
diffusion  operator  and  the  filtering  operator  on  the  same  grid.  The  results 
on  the  right  were  obtained  by  applying  the  filter  operator  on  a  grid  with 
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Figure  3:  True  solution  (left),  unstabilized  solution  (right) 


twice  the  resolution  of  the  solution  grid.  The  results  are  enhanced  by  filter¬ 
ing  on  the  finer  grid  because  the  goal  of  the  filter  is  to  resolve  errors  on  the 
subgrid  scale. 


Figure  4:  Smoothed  solution  on  one  grid  (left)  and  two  grids  (right) 


6  Temporal  Integration  for  Non-Newtonian  Flows 

We  have  done  a  significant  amount  of  work  has  been  done  to  improve  tem¬ 
poral  integration  for  non-Newtonian  flows.  The  modeling  equations  for 
non-Newtonian  flows  require  resolution  of  the  fluid’s  stress,  velocity,  and 
pressure.  The  conservation  equations  for  mass  and  momentum  represent  a 
”  Stokes-like”  system,  and  the  constitutive  equation  for  stress  is  a  hyperbolic 
partial  differential  equation.  Instead  of  solving  the  entire  coupled  system 
implicitly,  we  computed  using  sub-steps  of  a  given  time  step  At.  The  con¬ 
servation  equations  were  resolved  in  the  first  fraction  of  the  time  step,  using 
a  lagged  value  for  stress.  These  solutions  were  used  in  the  next  fraction  of 
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the  step  to  update  the  stress,  and  the  updated  value  for  stress  was  used 
to  compute  a  final  velocity  and  pressure  over  the  time  step.  This  strategy 
allows  one  to  apply  stabilization  techniques  over  a  fraction  of  a  time  step  in¬ 
stead  of  the  entire  time  step,  reducing  the  damping  effects  of  added  diffusion 
operators. 

We  began  the  research  by  focusing  on  the  convection-diffusion  equation. 
While  operator  splitting  techniques  have  frequently  been  applied  to  this 
equation,  our  strategy  differed  in  that  we  combined  features  of  both  ad¬ 
ditive  decomposition  algorithms  (such  as  alternating  direction  implicit  and 
implicit/explicit  schemes)  with  product  decomposition  algorithms.  The  ad¬ 
ditive  decomposition  algorithms  decompose  the  operator  and  treat  variables 
either  implicitly  or  explicitly  at  each  fraction  of  the  time  step.  Product 
decomposition  algorithms  advance  the  convective  part  of  the  operator  first 
and  use  this  to  advance  the  diffusive  part  of  the  operator  over  the  entire 
time  step. 

We  obtained  optimal  a-priori  estimates  for  our  scheme  for  both  the 
convection-diffusion  equation  and  the  equations  for  viscoelasticity.  We  showed 
the  validity  of  the  scheme  for  several  applications  problems,  including  a  mov¬ 
ing  front  problem,  a  rotating  cone  problem,  and  a  contraction  problem  from 
the  viscoelasticity  literature.  John  Chrispell  received  his  Ph.D.  in  2008,  un¬ 
der  my  supervision  and  that  of  Professor  Vince  Ervin.  Dr.  Chris  Kees  of 
the  US  Army  ERDC  served  on  his  Ph.D.  committee.  John  spent  summers 
at  ERDC  and  implemented  his  algorithm  in  PyADH,  an  ERDC  research 
hydraulics  code. 

The  work  on  this  project  was  jointly  funded  by  NSF  under  grant  No. 
DMS-0410792  and  this  ARO  grant.  It  is  summarized  in  [3,  4,  5]  and  [2], 

7  Computation  of  Navier-Stokes  Equations 

V.J.  Ervin  has  worked  on  extending  the  enhanced-physics  based  scheme 
for  three-dimensional  Navier-Stokes  equations  [1].  The  general  enchanced 
physics  based  scheme  (preserving  helicity)  for  the  NS  equations  introduced 
by  Rebholz  is  based  on  expressing  the  NS  equations  in  curl-form.  In  the 
curl-fornr  the  pressure  term  in  the  equations  is  the  Bernoulli  pressure  ( p  + 
1/2|k|2).  Due  to  the  behavior  of  |rt|2  in  or  near  the  boundary  layers,  the 
Bernoulli  pressure  term  can  contribute  significant  computational  instabil¬ 
ity  to  the  approximation.  In  [1].  Ervin,  Rebholz,  et.al.  investigated  and 
showed  the  positive  effect  of  grad-div  stabilization  on  the  enhanced  physics 
approximation  scheme. 
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Ervin  is  also  working  in  improvements  to  computational  solutions  of  the 
Navier-Stokes  equations  using  filtering.  The  scheme  mirrors  the  scheme  out¬ 
lined  in  [6] ,  in  that  the  filtering  operator  acts  as  a  post-processing  operator 
for  the  solution  to  the  discrete,  unstabilized  operator.  There  are  theoretical 
results  to  show  that  the  filtering  scheme  reduces  the  effects  of  high-frequency 
components  that  appear  in  the  computed  solution 

The  paper  associated  with  this  work  will  be  submitted  within  the  next 
six  months. 

8  Student  Involvement 

While  student  involvement  does  not  directly  connect  with  scientific  progress, 
I  wanted  to  add  this  section  to  discuss  the  impact  that  this  funding  has  had 
on  John  Chrispell,  who  received  a  Ph.D.  in  mathematical  sciences  while 
being  partially  funded  by  ARO.  John  was  able  to  spend  two  summers  work¬ 
ing  with  Dr.  Chris  Kees  at  the  US  Army  ERDC  in  Vicksburg,  MS.  Chris 
is  an  expert  in  numerical  methods  for  porous  media  flow,  and  he  is  well 
versed  in  a  variety  of  software  languages.  John  learned  of  physical  ap¬ 
plications  for  his  work  on  convection-diffusion  equations,  he  enhanced  his 
computational  skills,  and  he  implemented  his  operator-splitting  algorithm 
in  a  real-world  simulation  code.  He  worked  in  an  interactive,  laboratory 
environment  and  was  able  to  learn  from  engineers  whose  daily  work  requires 
extensive  knowledge  of  flow  problems,  resolution  of  computed  solutions  with 
underlying  physical  behavior,  computational  skills,  and  the  ability  to  effec¬ 
tively  communication  through  papers  and  presentations.  It  is  impossible  to 
underestimate  the  benefit  he  received  from  his  work  there.  Students  whose 
whole  tenure  is  spent  in  the  purely  academic  setting  can  never  truly  embrace 
the  concept  of  ”  applied  mathematics” .  I  am  extremely  grateful  to  ARO  for 
allowing  me  to  provide  this  experience  for  my  students. 
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